arXiv:1507.03220vl [astro-ph.HE] 12Jul2015 


Shocks in unmagnetized plasma with a shear flow: Stability and magnetic fleld 

generation 

M. E. Dieckmann,^ A. Bock,^ H. Ahmed,^ D. Doria,^ G. Sarri,^ A. Ynnerman,^ and M. Borghesi^ 

^Department of Science and Technology, Linkoping University, SE-60174 Norrkoping, Sweden 
^Centre for Plasma Physics (CPP), Queen’s University Belfast, BT7 INN, UK 

(Dated: July 14, 2015) 

A pair of curved shocks in a collisionless plasma is examined with a two-dimensional particle-in¬ 
cell (PIC) simulation. The shocks are created by the collision of two electron-ion clouds at a speed 
that exceeds everywhere the threshold speed for shock formation. A variation of the collision speed 
along the initially planar collision boundary, which is comparable to the ion acoustic speed, yields a 
curvature of the shock that increases with time. The spatially varying Mach number of the shocks 
results in a variation of the downstream density in the direction along the shock boundary. This 
variation is eventually equilibrated by the thermal diffusion of ions. The pair of shocks is stable 
for tens of inverse ion plasma frequencies. The angle between the mean flow velocity vector of the 
inflowing upstream plasma and the shock’s electrostatic field increases steadily during this time. 

The disalignment of both vectors gives rise to a rotational electron flow, which yields the growth of 
magnetic field patches that are coherent over tens of electron skin depths. 

PACS numbers: 52.35.Tc 52.35.Fp 52.65.Rr 


I. INTRODUCTION 

The collision of an ionized blast shell with an ambient 
plasma triggers the formation of shocks if the collision 
speed exceeds a threshold value. The critical speed de¬ 
pends on the plasma wave mode that is mediating the 
shock and on the importance of Coulomb collisions be¬ 
tween particles. If the mean frequency, with which the 
plasma particles collide, is well below all resonance fre¬ 
quencies of the plasma, then the effects of binary colli¬ 
sions are negligible and the shocks are mediated by elec¬ 
trostatic and electromagnetic helds. 

The plasma processes that sustain a collisionless shock 
and the structure of the associated electromagnetic helds 
vary strongly between the different plasma regimes. So¬ 
lar system shocks, like the Earth’s bow shock mm , are 
immersed in a plasma that is carrying a relatively strong 
magnetic held and they connect two plasmas that col¬ 
lide at a non-relativistic speed. Such shocks are usually 
mediated by magnetosonic waves. The collision speed 
between a supernova blast shell and the ISM at a late 
evolution phase is similar to the collision speed between 
the solar wind and the Earth’s bow shock. The ampli¬ 
tude of the magnetic held in the ISM, into which a su¬ 
pernova remnant (SNR) shock expands, is weaker by an 
order of magnitude than that in the solar wind plasma 
at the Earth’s orbit. Magnetosonic waves, which have a 
low held amplitude, may not be able to sustain perma¬ 
nently the shock because other instabilities develop faster 
and on a smaller spatial scale. Simulations have shown 
that drift instabilities and electrostatic turbulence can in 
some cases suppress the growth of a magnetosonic wave 
[3]. As we go to higher how speeds, the plasma shocks 
become magnetized by hlamentation instabilities HE]- 

We consider here shocks, which develop in an initially 
unmagnetized and collisionless plasma. Such shocks are 
frequently observed in the laboratory [HHS] and they 


might be representative for SNR shocks in their late evo¬ 
lution phase. An electrostatic shock in its most basic 
form is characterized by a potential difference, which 
is sustained self-consistently by the plasma. The shock 
connects the downstream region and an upstream region 
ahead of the shock. Electrons stream from the denser 
downstream region into the upstream region and create 
a charge imbalance between both regions. The denser 
downstream plasma goes on a positive potential relative 
to the upstream plasma. 

The upstream plasma streams towards the shock at a 
speed, which exceeds the ion acoustic speed. The up¬ 
stream ions are slowed down and compressed by the po¬ 
tential jump as they cross the shock. The potential jump 
reflects some of the incoming upstream ions, which then 
move back upstream. The remainder of the incoming 
ions enters the downstream region, which expands due 
to the accumulation of the inflowing ions. The shock is 
thus not stationary in the downstream frame of reference 
and moves upstream. 

An electrostatic shock is an ion phase space structure, 
which consists of inflowing upstream ions, reflected ions 
and ions that overcame the positive potential and accu¬ 
mulated downstream of the shock. The electrostatic field, 
which mediates the shock, will also accelerate some of the 
downstream ions into the upstream direction and act as 
a double layer. Many shocks in unmagnetized plasma are 
a combination of a double layer and of an electrostatic 
shock. Such hybrid structures |I4j are characterized by 
a unipolar electric field. In what follows we shall refer to 
these hybrid structures as plasma shocks to distinguish 
them from pure electrostatic shocks. 

The incoming upstream ions, which have been reflected 
by the electrostatic shock, and the downstream ions, 
which have been accelerated upstream by the double 
layer, form a beam that outruns the plasma shock. We 
shall refer to this beam as the shock-reflected ion beam. 
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The number density of this beam can be a significant 
fraction of that of the incoming upstream ions, which 
implies that the ion number density ahead of the plasma 
shock is well above that in the far upstream region. In 
what follows, we refer to this region as the foreshock. 

The counterstreaming nonrelativistic and unmagne¬ 
tized ion beams drive the electrostatic ion acoustic insta¬ 
bility in the foreshock dsnn]. The speed of the shock- 
reflected ions in the upstream frame of reference exceeds 
the ion acoustic speed. The ion acoustic instability can, 
however, only be destabilized if the beam velocity com¬ 
ponent along the wave vector is subsonic. The wave vec¬ 
tors of the unstable waves can thus not be parallel to the 
plasma flow velocity vector. Oblique electrostatic waves 
grow and the obliquity angle is such that the beam ve¬ 
locity component along the wave vector is comparable 
to the ion acoustic speed m Obliquely propagating ion 
acoustic waves grow in the foreshock region and modulate 
the incoming upstream ions. The plasma shock is either 
transformed into a shock with a broad transition layer [3] 
or it is destroyed by the inflowing turbulent plasma [Si- 

Previous simulation studies have addressed the evo¬ 
lution of (quasi-)planar plasma shocks. The planarity 
has been enforced in the PIC simulations by resolving 
only one spatial direction or by choosing initial conditions 
that are uniform along one direction. However, in par¬ 
ticular the plasma shocks in laboratory experiments are 
often nonplanar. This motivates our study of the forma¬ 
tion and evolution of curved shocks with a particle-in-cell 
(PIC) simulation. 

The shock curvature is introduced in our simulation 
through the following setup. The two electron-ion clouds 
collide at a boundary, which is orthogonal to the spatially 
uniform collision direction at the simulation’s start. The 
mean speed of the plasma along the collision direction 
varies as a function of the orthogonal direction and this 
velocity shear gives rise to a shock front that becomes in¬ 
creasingly curved in time. The amplitude of the velocity 
change is comparable to the ion acoustic speed. 

We find that the shock formation and its stability are 
not affected by this large velocity shear. The life-time of 
the plasma shocks is of the order of tens of inverse ion 
plasma frequencies. The shock transition layer is trans¬ 
formed after this time by the onset of ion acoustic tur¬ 
bulence in the foreshock. The transition from a sharp 
electron skin depth-scale structure into a broad transi¬ 
tion layer is also observed for planar shocks. The key 
difference between the structure of the curved shock and 
a planar shock is tied to the disalignment of the elec¬ 
tric field with the flow velocity vector of the incoming 
upstream plasma. The disalignment gives rise to a rota¬ 
tional component of the electron flow, which yields the 
growth of magnetic field patches. These patches are co¬ 
herent over tens of electron skin depths and their size is 
limited by the simulation box dimension. Their ampli¬ 
tude yields a ratio of the electron plasma frequency to 
the cyclotron frequency of about 100. 

Our paper is structured as follows. Section 2 discusses 


the PIC code and the initial conditions we use. Section 3 
examines the formation and the evolution of the pair of 
electrostatic shocks. Section 4 summarizes our findings. 

II. INITIAL CONDITIONS AND THE PIC 
METHOD 

Particle-in-cell (PIC) simulation codes [T5] solve the 
Vlasov-Maxwell system of equations via the method 
of characteristics m- The electromagnetic fields are 
evolved in time via Ampere’s law and Faraday’s law. 
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Moeo-^ = V X B - ^oJ, (1) 


Most codes fullfill the equations V • E = p/eg and V • 
B = 0 either as constraints or via correction steps. The 
plasma is approximated by an ensemble of computational 
particles (CPs), which correspond to volume elements of 
the phase space density distribution. The charge-to-mass 
ratio of the CPs equals that of the plasma particles they 
represent. Their momentum and position are updated 
with the relativistic Lorentz force equation 

^ = (E(xj) -k vj X B(xj)). (3) 

The CP with the index j of the species i with the posi¬ 
tion Xj and the relativistic momentum pj = rriiTj'Vj has 
the charge qi and the mass rrii. Its position is updated 
by dbcj/dt = Vj. The CPs and the electromagnetic fields 
E(x, t) and B(x,<) are connected as follows. The fields 
are interpolated to the particle position and update its 
momentum in time. The current density contribution of 
each CP is interpolated to the grid. The summation over 
all current density contributions yields the macroscopic 
current density J(x, t), which updates the electromag¬ 
netic fields via Ampere’s law. 

A pair of shocks is generated in the simulation by the 
collision of two spatially uniform plasma clouds of equal 
density. Each cloud consists of electrons with the charge 
—e (e: elementary charge), the mass me, the number 
density no and the temperature Te = I keV. We model 
Deuterium ions with the number density no, the mass 
mu and the temperature To = 200 eV. We motivate our 
choice for the initial conditions as follows. 

The collision of two plasma clouds in laboratory- or 
astrophysical settings usually involves two plasmas with 
densities that differ by orders of magnitude. The blast 
shell, which is ejected during a supernova, is composed of 
stellar material. Its density is thus initially much higher 
than that of the ambient plasma, which is the stellar 
wind the star emanated prior to the supernova. A laser¬ 
generated blast shell is also much denser than the am¬ 
bient medium, which is the residual gas that has been 
ionized by secondary x-ray radiation from the target. 
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The blast shells are practically unaffected by the ambi¬ 
ent medium during their initial expansion phase and they 
expand in the form of rarefaction waves. The front of the 
blast shell becomes faster and thinner in time. Once the 
expansion speed of the blast shell becomes supersonic 
and once the ram pressure of the expanding blast shell 
becomes comparable to the thermal pressure of the ambi¬ 
ent medium, the front can be confined and shocks form. 
Experimental observations m and PIC simulations [20] 
suggest that shocks form when the densities of the collid¬ 
ing clouds become comparable. Selecting equal densities 
for both colliding plasma clouds should thus be a valid 
initial condition. 

The electron and ion temperatures we use are typical 
for plasmas, which are created when an ultrashort laser 
pulse ablates a solid target and if time scales are con¬ 
sidered, which are short compared to the time it takes 
to establish a thermal equilibrium via collisions between 
electrons and ions. The temperature of the electrons in 
the ambient plasma, which has been ionized by secondary 
x-ray emissions from the laser-ablated target, is conpara- 
ble to 1 keV. The temperatures of the ions of the ambient 
medium and of the blast shell are usually well below that 
of the electrons m- 

Deuterium ions have the same charge-to-mass ratio as 
the fully ionized light atoms, which we usually encounter 
in laser-plasma experiments, and they have the largest 
thermal velocity spread for a given temperature. Hence 
they will provide the strongest ion Landau damping. If 
we observe the growth of ion acoustic waves for counter¬ 
streaming beams of Deuterium, then the same will be 
true for equally hot plasmas, which consist of heavier 
ions with the same charge-to-mass ratio as Deuterium. 

The plasma frequencies of the electrons and of the ions 
are Wp,e = (noe^/eoWe)^^^ and ujp^i = 
respectively. The ion acoustic speed in this plasma is 
Cs = {ickBiTf,+Tu)/mD) ' assuming that the adia¬ 
batic constant 7 c = 5/3 is the same for both species. 
The ion acoustic speed is Cs = 3.1 • 10® m/s. 

The simulation domain has the dimensions x Ly = 
176As X 26.2As where A^ = c/wp^e is the electron skin 
depth. The boundary conditions are periodic along y 
and open along x. Choosing periodic boundaries along 
y implies that our simulation evolves in time a periodic 
chain of blast shells rather than a solitary one. The sim¬ 
ulation domain is split up into two equal parts along x. 
The blast shell occupies the interval < x < 0 

and —Lp/2 < y < Ly/2. The second cloud, which we 
refer to as the ambient plasma, occupies the interval 
0 < a; < La;/2 and —Ly/2 < y < Lj//2. 

The initial mean speeds of the electrons and ions of 
the ambient plasma vanish. The mean speed of the blast 
shell’s electrons equals that of the ions and is denoted 
here as v = (■!;2,,0,0). The value of Vx varies piecewise 
linearly along the y-axis. The largest value ofvx = 1.15 x 
10® m/s or 3.7cs is reached at the position y = 0. The 
speed decreases linearly in both y-directions. It reaches 
its minimum of 8.7 x 10® m/s or 2.8cs at the boundaries 


at y = —Lp/2 and y = Ly/2. 

The simulation box is resolved by 4000 grid cells along 
the x-direction and by 600 cells along the y-direction. 
The quadratic side length of each cell is Ax = 0.044As. 
Each plasma species is resolved by 200 CPs per cell. We 
evolve the system for = 153 using 3.2 x 10® 

time steps A(. In what follows, we normalize time to 
l/wpy, space to As and speed to the electron thermal 
speed Vth,e = {ksTe/me)^^'^■ The electric field is normal¬ 
ized to mf^LOp^ec/e and the magnetic one to meWp^e/e. 

III. THE SIMULATION RESULTS 

In what follows, we shall present and discuss the spa¬ 
tial distributions of the ion density, of the amplitude of 
the flow-aligned electric field component Ex{x, y) and the 
out-of-plane magnetic field distribution Bz{x,y) at the 
times t = 6, 13.8, 19.4, 33.8, 70 and 153. 

Figure]^ shows the ion density, the electric Ex compo¬ 
nent and the magnetic component close to the initial 
contact boundary a; = 0 at several times. A band with an 
increased ion density is visible in Fig. Ha). The ions of 
both clouds interpenetrate and their cumulative density 
exceeds Uq- The peak density of the ions is reached at 
cc « 0.7 and it exceeds 2no for all values of y. This ion 
cloud overlap layer is broadest at y « 0. At this time the 
ions of both clouds move independently and the width of 
the layer is proportional to the speed at which the clouds 
collide. The ion density is not constant close to its max¬ 
imum value in Fig. m- Hence a downstream region, 
which is characterized by a spatially uniform plasma dis¬ 
tribution along X that separates the forward and reverse 
shocks, has not formed at this time. 

A strong bipolar electric field pulse is visible in Fig. 
m- It is sustained by the space charge, which results 
from the electrons that escaped from the ion cloud over¬ 
lap layer. The polarity of the electric field is such that the 
ion cloud overlap layer is on a positive potential relative 
to the ambient- and blast shell plasmas. The slow-down 
of the inflowing ions by this potential is responsible for 
the increase of the ion density beyond 2no. 

Figure H c) reveals magnetic oscillations within the ion 
cloud overlap layer. The strongest oscillations are located 
at y « 0 and x « 0.7, where the plasma collides at the 
highest speed. The growth of B^ within the ion cloud 
overlap layer can be attributed partially to the instabil¬ 
ity observed by [22] for shocks and by |23| and |21| for 
rarefaction waves. The electrons are accelerated along 
the electric field of the ion cloud overlap layer. The re¬ 
sulting directional anisotropy in the velocity distribution 
triggers the Weibel instability [211 |2S| . 

The Weibel instability in the form discussed in Ref. 
[26] can, however, not explain why the magnetic field 
oscillations peak at y « 0 and x « 0.5. We attribute this 
to the geometric effect that is outlined in Fig. 

We consider first the boundaries between the over¬ 
lap layer and the upstream plasma that have a constant 
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FIG. 1: The ion density distributions expressed in units of no are shown in the first column (from left to right). The distributions 
of Ex{x,y) are shown in the second column. The electric amplitudes have been multiplied by a factor 100. The third column 
shows the distributions of B^{x,y) and the right column shows slices of the ion density distributions along y — 0 (dashed blue) 
and along y = 13.1 (black). The upper row corresponds to the time t — 6, the second row to t — 13.8, the third one to t = 19.4 
and the bottom row to t = 33.8. 


slope. Electrons that stream into the overlap layer are 
accelerated by the ambipolar electric held that ensheaths 
the overlap layer. The held is aligned with the boundary 
normal. The electrons are accelerated and the ions are 
decelerated when they enter the overlap layer, while their 
lateral velocity components remain unchanged. The ve¬ 
locity vectors of the inhowing electrons and ions are thus 
rotated into opposite directions. Their current contribu¬ 
tions do no longer cancel each other out and a net current 
develops within the overlap layer. 

An even stronger net current develops close to the 
cusps due to the changing direction of the normal. Elec¬ 
trons that cross the overlap layer at a concave cusp are 
scattered. Electrons that cross it at a convex cusp are fo¬ 


cused. The electron current is no longer balanced along 
the vertical direction in the center of Eig. A net cur¬ 
rent flows from the concave to the convex cusp of the 
overlap layer and a magnetic field will grow, which is at 
least initially confined to the overlap layer. The grow¬ 
ing spatially localized magnetic field will induce a ring 
current within the overlap layer. The direction of the 
magnetic field in Eig. [^matches that in the simulation if 
we take into account that the z-axis points into the plot’s 
plane in Fig. [^c). 

Magnetic held oscillations are present on both sides of 
the ion cloud overlap layer in Fig. [^c). Their wave¬ 
length is 2tt/L y along y. This modulation extends up 
to the boundary at —L^I2 and it does not oscillate 
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FIG. 2: The field and fiow diagram within the overlap layer 
of both plasma clouds. The normals of the boundary between 
the overlap layer and the upstream plasma are denoted by N 
and are parallel to the vector of the ambipolar electric field. 
Electrons are denoted by e and their velocity vectors are the 
dashed arrows. Deuterium ions are denoted by d and their 
velocity vectors are solid. Net currents are denoted by solid 
vectors and the symbol J. The circulating current gives rise 
to a magnetic field (Symbol B). 


along the x-direction (not shown). Such an oscillation 
can not be driven by an electron current that emanates 
from the overlap layer. Electrons that move at the ther¬ 
mal speed can only traverse the distance ~ 15As dur¬ 
ing t = 6. Hence the magnetic field must be driven 
by spatial and temporal variations of the electric field 
and propagate in the light mode. Faraday’s law gives 
■§f.Bz = -^Ex — -§^Ey. The observation of a magnetic 
field Bz ^ 0 indicates that the condition for a two- 
dimensional electrostatic potential dE^/dy = dEJdx 
is not fulfilled by the bipolar structure in Fig. Qb). 
The noise and the Weibel instability, which triggers the 
growth of patchy magnetic fields in the ion cloud overlap 
layer of the ions, introduce a weak magnetic component 
of the shock to start with. The magnetic fields, which 
grow upstream of both shocks, remain weak and we shall 
not discuss them further. 

Figures l^e-h) show the ion- and electromagnetic field 
distributions sampled at t = 13.8. The ion cloud over¬ 
lap layer in Fig. [^e) is now broader and less dense at 
the boundaries y = ±Ly/2 than in the center, which is 
demonstrated quantitatively by Fig. The ion den¬ 

sity in Fig. downstream of the slow shocks evidences 
a flat density profile in the interval 0.7 < x < 1.8, which 
is typical for a downstream region. 

The larger plasma compression by fast shocks com¬ 
pared to that of slow shocks implies that the density ra¬ 
tio between the downstream plasma and the upstream 
plasma is larger for the shocks close to y « 0 in Fig. 
[^e). Fast shocks like the ones at y « 0 reflect most of 
the incoming ions and only a minor fraction enters the 
downstream region [izi m\- Consequently, the down¬ 


stream region of the fast shocks expands slowly. More of 
the incoming upstream ions can traverse the slow shocks 
close to y = ±Ly/2 and the downstream region behind 
them accumulates more ions. It expands faster. The 
ion density in Fig. is decreasing rapidly to a value 
« l.Sng at y = ±Ty/2 as x is decreased below 0.7 or 
increased above 1.8. The density converges to « Uq at 
the boundaries of the displayed interval. The density en¬ 
hancements between — 1 < a; < 0.5 and 2 < a; < 4 are 
caused by the shock-reflection of ions and by ions that 
propagated through the ion cloud overlap layer before 
both shocks formed. 

The electric field outlines the location of the forward 
and reverse shocks in Fig. Sf). The unipolar electric 
field pulses, which mediate the shocks, are closest at 
?/ « 0. Their separation increases along x as we go to the 
boundaries at y = ±Lyj2. The thickness of the unipo¬ 
lar electric field peaks along x is less than in Fig. Bb). 
This change of the thickness of the pulse evidences the 
transformation of the ion cloud overlap layer into a down¬ 
stream region that is enwrapped by forward and reverse 
shocks. 

The magnetic field amplitude modulus in Fig. [^g) 
has quadrupled at y « 0 compared to that at the ear¬ 
lier time. Additional field structures have emerged close 
to the boundaries at y = ±Ly/2 and at y = ±Ly/A. 
The magnetic field is strongest inside of the downstream 
region. It is, however, not confined by it like in the sim¬ 
ulation in Ref. [H]. The strongest magnetic fields are 
observed in the intervals along y, where the cloud colli¬ 
sion speed has extrema and where the overlap layer has 
cusps. The normalization of the fields implies that the 
peak amplitude of the magnetic field yields a ratio of 
the electron cyclotron frequency to the electron plasma 
frequency of about 6 x 10“^. 

The ion distribution has changed qualitatively at the 
time t = 19.4, which is evidenced by the Figs. Bi-l)- 
The ion density in the downstream region is comparable 
to that at the earlier time t = 13.8, but the ion density 
is now also a function of x. The density of the down¬ 
stream ions in Fig. [^i) and in Fig. Bi) is highest at 
the concave shocks and lowest at the convex shocks. The 
ion density at y = 0 and a; « 2 is about 3no at the con¬ 
cave reverse shock and it decreases to a value of 2.7no 
at x « 2.5. The electric field in Fig. Bj) reacts to it 
because a larger change of the ion density across the con¬ 
cave shock yields a stronger ambipolar electric field. The 
electric field modulus at the concave shock at y = 0 and 
X « 1.8 exceeds that of the convex shock at x ~ 2.5. The 
magnetic B^ component in Fig. m component shows a 
cellular structure and peak amplitudes of « 8 x 10“^ are 
reached close to y « 0. 

Figures [^m-p) show the ion density distribution and 
the electromagnetic field distributions at the time t = 
33.8. The ion density distribution in Fig. B™) the 
two density slices shown in Fig. Bp) demonstrate that 
the density of the downstream plasma immediately be¬ 
hind a shock still depends on whether it is concave or 
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FIG. 3: The ion phase space density distribution fi{x,y,Vx) 
at the time t = 33.8. The blast shell ions (blue) are located 
mainly at high positive values of Vx while the ambient ions 
(red) are located mainly at low values of Vx- The x-axis range 
is [-2.75,11.0], the y-axis range is [-13.1,13.1] and the U 2 :-axis 
range expressed in units of Cs is [-1.9,5.8]. 


convex. The ion density decreases with increasing x close 
to y = 0, while the opposite is true close to y = ±Ly/2. 

The ion density distribution in the foreshock region of 
each shock has changed from a diffuse distribution in Fig. 
[^i) to one that shows a pronounced peak that is located 
about 2As ahead of each shock. These structures yield 
a thin band in Fig. [^n) with an amplitude modulus of 
about 5 X 10“^. The electric field of such a structure and 
the electric field of the nearest shock have the opposite 
polarization. 

The cellular magnetic field structures from Fig. m 
have merged to a large magnetic field distribution in Fig. 
[^o), which extends far into the upstream regions of both 
shocks. The magnetic field amplitude peaks in the down¬ 
stream region close to y « 0. A weaker similar distribu¬ 
tion exists close to the boundary at y = ±Lyl2. The 
shape of the magnetic field structure suggests that the 
associated current has its source in the kink in the over¬ 
lap layer at y = 0. We expect that the net current at 
this kink is higher than that at y = Ty/2, because the 
plasma flow speed is highest at y = 0 and because the 
ion density and, hence, the potential of the overlap layer 
peak at y = 0. 

The shock structure and the source of the ion density 
peaks ahead of the main shocks in Fig. [^m) is revealed 
by the ion phase space density distribution fi{x,y,Vx), 
which is displayed in Fig. This distribution shows 
several distinct features. The incoming blast shell ions 
are located at low values of x and at high positive speeds 
Vx- The blast shell ions have their largest mean speed at 
y « 0 and their speed decreases linearly as we go towards 
the periodic boundaries in the y-direction. The blast 
shell propagates to increasing values of x. The ambient 
ions are located at the right (large values of x) and their 
mean speed is zero. The ions of both clouds merge to a 
structure, which is characterized by a large spread along 
the Ujj-axis. This is the downstream ion population and 
it is bounded by the forward and reverse shocks along 
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FIG. 4: The ion phase space density distribution fi{x,y,Vx) 
at the time t = 33.8 and along the slice y = 0. Positions are 
normalized to As and the velocity axis is normalized to Cs. 
Panel (a) shows the total ion density. Panel (b) shows the 
blast shell ions and panel (c) the ambient ions. 


both ^-directions. The velocity change between the blast 
shell plasma and the downstream region and between the 
downstream region and the ambient plasma are caused 
by the ion acceleration by the shock’s electric held. 

Each plasma shock gives rise to a shock-reflected ion 
beam. This beam is composed of the ions that were re¬ 
flected by the electrostatic shock and of the ions that 
were accelerated by the double layer as they moved from 
the downstream region into the upstream region. Let us 
consider the shock-reflected ion beam, which is located 
to the left and at low speeds. This beam reveals two dis¬ 
tinct regions. The part to the left is uniform along the 
y-direction and it contains only ions from the ambient 
plasma. Its mean velocity and its density decreases as 
we go to decreasing values of x. Such a phase space pro¬ 
file is that of a rarefaction wave [23] . The shock-reflected 
ion beam close to the shock is no longer spatially uniform 
along y. The boundary between these two regions follows 
the shock profile and it separates the rarefaction wave, 
which contains only ions from the ambient plasma, from 
the shock-reflected ion beam that contains ions from the 
blast shell plasma and the ambient plasma. The shock- 
reflected ion beam at large values of x and at positive Vx 
also shows a clear subdivision into two domains, which 
are separated by a boundary. The mean velocity of the 
ions changes across this boundary and the ion acceler¬ 
ation is accomplished by the electric field pulse seen in 
Fig. [^n) at large x. 

A slice of the ion phase space distribution shown in Fig. 
along X and for y = 0 is displayed in Fig. Figure Qa) 
shows the cumulative ion distribution. The blast shell 
ions are found in the interval cc < 3 and Vx/Cs ^ 4: and the 
ambient ions at a: > 5 and Vxjcs ~ 0. The downstream 
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FIG. 5: The ion density distribution expressed in units of 
no is shown in panel (a). Panel (b) shows the distribution 
of Ex{x,y). The ion density distributions along two slices 
y = 0 (dashed black) and y — 13.1 (black) are shown in panel 
(c). Panel (d) shows the distribution of B^{x,y). The time is 
t = 70. 


region is confined to 3.5 < x < 4.5. The blast shell 
ions in Fig. |^b) are partially reflected at a; « 3.5 and 
some of them enter the downstream region. This phase 
space structure is an electrostatic shock. Some of the 
ions have crossed the downstream region and they have 
been accelerated by the double layer at a; « 4.5. The 
beam of accelerated ions has an almost constant speed 
up to a; « 6, where the beam speed decreases. The beam 
speed grows linearly and the beam density decreases as 
we go from x « 6.5 to x « 10.5. The latter phase space 
structure is a rarefaction wave. The sudden decrease of 
the mean speed of the ion beam at a; « 6 in Fig. ib) 
implies that ions at lower x catch up and collide with 
ions at higher values of x. 

Figure shows the ion density distribution and the 
electromagnetic field distributions at the time t = 70. 
The downstream region close to y « 0 is no longer 
bounded by smooth shocks. An ion density cusp be¬ 
tween the shock and the blast shell has developed at 
a; « 7 and y « 0, while the shock is convex on the 
other side of the downstream region and surrounded by 
cusps at y = ±2. The fastest shocks show a more com¬ 
plex distribution than their slower counterparts close to 
y « ±Ly/2. The ion density distribution along x in Fig. 
[^c) of the latter is qualitatively similar to that at previ¬ 
ous times; the ion density at the concave shocks exceeds 
the one at the convex shocks. The amplitude modulus 
of the electric field peaks, which mediate both shocks, 
has decreased from a value « 2 x 10“^ at t = 33.8 to 
an amplitude modulus of about 1.5 x 10“^ at t = 70. 
The extent of the magnetic field patches in Fig. id) is 
now of the order of 10 Ag. Any further expansion of the 
magnetic field patches along y is impeded by the periodic 



FIG. 6: The ion phase space density distribution fi{x,y,Vx) 
at the time t = 70. The blast shell ions (blue) are located 
mainly at high positive values of while the ambient ions 
(red) are located mainly at low values of Vx- The x-axis range 
is [-5.75,12.0], the y-axis range is [-13.1,13.1] and the U 2 ,-axis 
range expressed in units of Ca is [-1.9,5.8] (multimedia view). 


boundary conditions along this direction. 

The ion phase space density distribution at t = 70 in 
Fig. § looks qualitatively similar to that at the earlier 
time t = 33.8. The ions have started to mix in the down¬ 
stream region enclosed by both shocks. This mixing is 
accomplished by the ion phase space vortices, which we 
can observe close to the interface between the ambient 
and blast shell ions. The shock-reflected ions have prop¬ 
agated farther away from the shocks, but there is still a 
clear subdivision into ions, which were accelerated by the 
double layers before the shocks formed, and ion beams 
that were accelerated after the shock formation. The 
latter consist of ions of the blast shell plasma and of the 
ambient ions. The shock-reflected ambient ions to the 
right of the figure have started to overtake some of the 
blast shell ions. 

Figure shows the ion density distribution at t = 153. 
The downstream region in Fig. [^has a density that does 
no longer vary as a function of y. The density in the 
entire overlap layer is about that observed earlier close 
to the fastest shocks. Thermal diffusion is one way to 
equilibrate the ion density in the downstream region. The 
ion diffusion length can be estimated by multiplying the 
simulation time t = 153 with the initial thermal speed 
« 10^ m/s of the ions. Ions moving at this speed can 
cross the distance 5m ~ Ly/%. The shock-heated ions 
in the downstream region have speeds well in excess of 
the initial thermal speed and they can cross the distance 
from the high density region at y = 0 to the low density 
region at y = ±Lj,/2 during the simulation time. The 
equilibration of the ion density in the downstream region 
can thus be accomplished by the thermal diffusion of ions. 

Figure shows the distribution of the electric Aj, com¬ 
ponent at t = 153. Figure [^reveals that the shock tran¬ 
sition layer has changed from being a narrow unipolar 
electric field pulse observed at t = 70 to a broad layer of 
electrostatic waves. Such a turbulent layer is driven by 
ion acoustic waves, which reach electric field amplitudes 















FIG. 7: The ion density distribution expressed in units of no 
at the time t = 153 (multimedia view). 
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FIG. 8: The spatial distribution of the electric i5a;-component 
at t = 153 (multimedia view). 


that are about 50% of that of the electrostatic shock in 
Fig. [^b) . Such a turbulence layer is capable of thermal- 
izing the incoming ions in all directions, since the ions 
are exposed to a series of strong electric field pulses with 
an almost random polarization. The broader shock tran¬ 
sition layer also results in a lower ion density gradient in 
Fig.0 A decrease in the ion density gradient brings with 
it a low amplitude of the ambipolar electric field. That 
is the reason for why we can no longer see the shock’s 
unipolar electric field even though the potential differ¬ 
ence between the downstream region and the upstream 
region has not changed compared to that at t = 70. 

Figure shows the distribution of the magnetic 
component at t = 153. The magnetic field patches in Fig. 
I^have expanded along x compared to those at t = 70 and 
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FIG. 9: The spatial distribution of the magnetic B^- 
component at t = 153 (multimedia view). 



FIG. 10: The ion phase space density distribution fi{x,y,Vx) 
at the time t = 141. The blast shell ions (blue) are located 
mainly at high positive values of Vx while the ambient ions 
(red) are located mainly at low values of Vx- The white band 
shows phase space intervals that are occupied by ions from 
the blast shell plasma and from the ambient plasma. The x- 
axis range is [-8.8,22], the y-axis range is [-13.1,13.1] and the 
Ux-axis range expressed in units of Cs is [-1.9,5.8]. 


their width along this direction is about dOA^. Their ex¬ 
pansion along the y-direction was already limited by the 
simulation box size at t = 70 and hence the patches could 
not expand further along this direction. The coherence 
scale of the magnetic field patches, their expansion far 
upstream of the shock and their close correlation with 
the cusps of the overlap layer exclude the Weibel insta¬ 
bility as the cause. The magnetic fields driven by the 
Weibel instability oscillate in space and their wavelength 
is comparable to an electron skin depth. 

The ion phase space density distribution at t = 141 is 
shown in Fig. This distribution demonstrates that a 
downstream region still exists. The shocks, which enclose 
this region, are more diffuse than at the earlier times. 
This is a consequence of the ion acoustic turbulence, 
which is now mediating the shocks. The interaction of 
the blast shell ions and the ambient ions with the turbu¬ 
lent wave fields mixes both populations in phase space. 
The phase space interval, in which both ion species have 
mixed, is indicated with the white band in Fig. |10| The 
white interval corresponds to voxels, in which we find 
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ions from both plasma clouds. 


IV. SUMMARY 

We have examined the formation and evolution of a 
pair of shocks in an initially unmagnetized plasma. Two 
plasma clouds, which consisted of electrons and ions, col¬ 
lided at a boundary, which was initially planar. The elec¬ 
trons and ions of both clouds had the same density. The 
electrons of both clouds had the same temperature. The 
electron temperature exceeded that of the ions by a fac¬ 
tor of 5. The electrons and ions of each cloud had the 
same mean speed at any position and the plasma was 
initially free of any net charge and current. The collision 
speed normal to the initial collision boundary varied as 
a function of the position along the boundary. The colli¬ 
sion speed was highest in the center of the simulation box 
and decreased linearly with a decreasing distance to the 
periodic boundary, on which the collision speed reached 
its minimum value. 

We have obtained the following results. A pair of 
shocks can form in a plasma with a large velocity shear. 
These hybrid structures that consist of an electrostatic 
shock and of a double layer have a lifetime that is compa¬ 
rable to that of planar shocks. The shock normal is ini¬ 
tially anti-parallel to the velocity vector of the incoming 
upstream ions. The incoming upstream ions are slowed 
down and compressed along the shock normal direction 
and no particle acceleration takes place in the plane that 
is orthogonal to the shock normal. The plasma dynamics 
involves only the position and velocity along the shock 
normal and the shock dynamics is one-dimensional. Con¬ 
sequently the variation of the shock speed in the down¬ 
stream frame of reference and the plasma compression 
with the collision speed resembles that in the parametric 
study of one-dimensional plasma shocks of [55]. 

The shock dynamics becomes two-dimensional after 
about 70 inverse ion plasma frequencies. The ion density 
differences in the downstream region are equilibrated by 
thermal diffusion. The density of the downstream region 
equilibrates at the previously highest value, which was 
reached behind the fastest plasma shock. The transition 
layer of the plasma shock is transformed from a sharp 
unipolar electric field pulse into a broad layer of electro¬ 
static turbulence. The latter converts the directed flow 
energy of the upstream ions into heat along all three ve¬ 
locity directions, which leads to a full thermalization of 
the inflowing upstream ions by the shock. 

The formation of the plasma shock and the associated 
increase of the ion compression yields a sudden increase 
of the positive potential of the double layer. The ions 
that cross the double layer towards the upstream direc¬ 
tion are accelerated to a higher speed after the shock has 
formed and they catch up with the ions that crossed the 
double layer at an earlier time. The magnitude of the 
velocity change is here sufficient to trigger the formation 
of a shock, which is located ahead of the main shock and 


has a life-time of the order of a few inverse ion plasma 
frequencies. 

The formation of a secondary shock ahead of the pri¬ 
mary one has been observed by [5^ . The secondary shock 
has been attributed the heat wave, which outran the ra¬ 
diative shock in this experiment. We can compare the 
life-time of this secondary shock to that we have observed 
in our simulation. The residual gas in Ref. |55| consisted 
of a mixture of Xenon and Nitrogen gas with a mass 
density of 3.6 x 10“^gcm“^. Let us consider the case 
study in Ref. [53], where the residual gas consists en¬ 
tirely of nitrogen and we furthermore assume that the 
nitrogen is fully ionized. We obtain an ion plasma fre¬ 
quency Wi « 3 X Our simulation would cover a 

time scale of about 50 ps, which is more than three or¬ 
ders of magnitude shorter than the life-time of the second 
shock in Ref. ]59]. A lower ionization state of the nitro¬ 
gen and the presence of neutral nitrogen would reduce 
the ion plasma frequency and extend the life-time of the 
secondary shock. A mix of nitrogen with different ioniza¬ 
tion states may also affect this life-time. It is, however, 
unlikely that the life-time of the secondary shock could 
be extended by a factor of 1000. The secondary shock, 
which has been observed in Ref. [53] , can thus not be ex¬ 
plained in terms of the subshock we have observed here. 

We observed the growth of large magnetic field patches 
in our simulation. Initially the magnetic field growth was 
limited to the ion cloud overlap layer. Magnetic fields 
can be generated via the Weibel instability [25] in spa¬ 
tially localized ion density accumulations such as shocks 
[55] and rarefaction waves [55] ]53] . The magnetic field 
structures observed at later times expanded into the up¬ 
stream region and they were not showing spatial oscilla¬ 
tions on an electron skin depth-scale, which are typical 
for the magnetic fields driven by the Weibel instability. 
The magnetic fields were coherent over tens of electron 
skin depths and the area they covered was limited by 
the dimensions of the simulation box and by the simula¬ 
tion time. The large-scale magnetic fields started to grow 
when the shock normal was no longer aligned with the 
plasma flow velocity vector. We have attributed the large 
scale magnetic field to currents, which initially develop 
close to the cusp in the overlap layer. The simulation 
shows that the magnetic field eventually diffuses out of 
the overlap layer. 

Experimental observations indicate that some SNR 
shocks are immersed in magnetic fields with amplitudes 
that exceed by far the values one would expect from the 
shock compression of the magnetic held of the interstellar 
medium |30| . Cosmic rays can magnetize the interstellar 
medium on large spatial scales [5T]. We scale the growth 
time of the magnetic held and the size of the magnetic 
patches to the plasma parameters found close to SNR 
shocks in order to determine if the corrugation of plasma 
shocks could be important for the magnetic held gener¬ 
ation at SNR shocks. We take the reference value 10 
cm“^ for the ion number density close to an SNR shock. 
Our simulation duration would correspond to » 4 x 10“^ 
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s. The spatial size « 40As of the magnetic field patches 
would correspond to about 50 km and their field ampli¬ 
tude would be of the order of 10 nT. The values for the 
growth time and the size of the magnetic field patches are 
microscopic compared to the size and the evolution time 
of an SNR shock. However, the magnetic field amplitude 
generated in our simulation exceeds that of the inter¬ 
stellar magnetic field by an order of magnitude. A cor¬ 
rugated shock front could thus generate magnetic fields 
ahead of the shock which are significantly stronger than 
those of the interstellar medium and it could compress 
these as it propagates across them. 

The periodic boundary conditions along the y-direction 
have limited the lateral expansion of the magnetic field 
patch at late times. An electron, which moves at the ther¬ 
mal speed, could cross the simulation box several times 
along the y-direction during the simulation time. Numer¬ 
ical artifacts, which are caused by the wrap-around of 
electrons, can usually be neglected because the electrons 
are scattered on the way by the electrostatic simulation 


noise. 

The key findings of this paper should however not be 
affected by the periodic boundary conditions. The net 
current that drives the magnetic field is generated in a 
small spatial interval close to the cusps that is far from 
the boundaries. The simulation box geometry will af¬ 
fect the shape of the generated magnetic field, but not 
its generation mechanism. The stability of the shocks is 
also not affected by the boundary conditions because the 
thermal speed of the ions is not high enough to let them 
cross the simulation box during the simulation time. 

Acknowledgements: The simulations were per¬ 
formed on resources provided by the Swedish National 
Infrastructure for Computing (SNIC) at HPC2N (Umea). 
GS acknowledges the EPSRC grant EP/L013975/1. 
The EPOCH code used in this research was devel¬ 
oped under UK Engineering and Physics Sciences Re¬ 
search Council grants EP/G054940/1, EP/G055165/1 
and EP/G056803/1. 


[1] S. D. Bale, M. A. Balikhin, T. S. Horbury, V. V. Kras- 
noselskikh, H. Kucharek, E. Mobius, S. N. Walker, A. 
Balogh, D. Burgess, B. Lembege, E. A. Lucek, M. Sc- 
holer, S. J. Schwartz, and M. F. Thomsen, Space Sci. 
Rev. 118 161 (2005). 

[2] D. Burgess, E. A. Lucek, M. Scholer, S. D. Bale, M. 
A. Balikhin, A. Balogh, T. S. Horbury, V. V. Kras- 
noselskikh, H. Kucharek, B. Lembege, E. Mobius, S. J. 
Schwarz, M. F. Thomsen, and S. N. Walker, Space Sci. 
Rev. 118 205 (2005). 

[3] M. E. Dieckmann, G. Sarri, D. Doria, H. Ahmed, and M. 
Borghesi, New J. Phys. 16 073001 (2014). 

[4] Y. Kazimura, F. Califano, J. 1. Sakai, T. Neubert, F. 
Pegoraro, and S. Bulanov, J. Phys. Soc. Jpn. 67 1079 
(1998). 

[5] J. T. Frederiksen, C. B. Hededal, T. Haugbolle, and A. 
Nordlund, Astrophys. J. 608 L13 (2004). 

[6] A. Spitkovsky, Astrophys. J. 673 L39 (2008). 

[7] A. Stockem, F. Fiuza, A. Bret, R. A. Fonseca, and L. O. 
Silva, Sci. Rep. 4 3934 (2014). 

181 D. W. Koopman, and D. A. Tidman, Phys. Rev. Lett. 
18 533 (1967). 

[9] S. O. Dean, E. A. McLean, J. A. Stamper, and H. R. 
Griem, Phys. Rev. Lett. 27 487 (1971). 

[10] A. R. Bell, P. Choi, A. E. Dangor, O. Willi, and D. A. 
Bassett, Phys. Rev. A 38 1363 (1988). 

[11] L. Romagnani, S. V. Bulanov, M. Borghesi, P. Audebert, 
J. C. Gauthier, K. Lowenbriick, A. J. MacKinnon, P. 
Patel, G. Pretzler, T. Toncian, and O. Willi, Phys. Rev. 
Lett. 101 025004 (2008). 

[12] G. Gregori et al, Nat. 481 480 (2012). 

[13] H. Ahmed, M. E. Dieckmann, L. Romagnani, D. Do¬ 
ria, G. Sarri, M. Cerchez, E. lanni, 1. Kourakis, A. L. 
Giesecke, M. Notley, R. Prasad, K. Quinn, O. Willi, and 
M. Borghesi, Phys. Rev. Lett. 110 205001 (2013). 

[14] N. Hershkowitz, J. Geophys. Res. 86 3307 (1981). 

[15] H. Karimabadi, N. Omidi, and K. B. Quest, Geophys. 


Res. Lett. 18 1813 (1991). 

[16] T. N. Kato, and H. Takabe, Phys. Plasmas 17 032114 

( 2010 ). 

[17] D. W. Forslund, and C. R. Shonk, Phys. Rev. Lett. 25 
281 (1970). 

[18] J. M. Dawson, Rev. Mod. Phys. 55 403 (1983). 

[19] T. H. Dupree, Phys. Fluids 6 1714 (1963). 

[20] G. Sarri, G. C. Murphy, M. E. Dieckmann, A. Bret, K. 
Quinn, 1. Kourakis, M. Borghesi, L. O. C. Drury, and A. 
Ynnerman, New J. Phys. 13 073023 (2011). 

[21] K. Eidmann, J. Meyer-ter-Vehn, and T. Schlegel, Phys. 
Rev. E 62 1202 (2000). 

[22] A. Stockem, T. Grismayer, R. A. Fonseca, and L. O. 
Silva, Phys. Rev. Lett. 113 105002 (2014). 

[23] C. Thaury, P. Mora, A. Heron, and J. C. Adam, Phys. 
Rev. E 82 016408 (2010). 

[24] K. Quinn, L. Romagnani, B. Ramakrishna, G. Sarri, M. 
E. Dieckmann, P. A. Wilson, J. Fuchs, L. Lancia, A. 
Pipahl, T. Toncian, O. Willi, R. J. Clarke, M. Notley, A. 
Macchi, and M. Borghesi, Phys. Rev. Lett. 108 135001 
( 2012 ). 

[25] E. S. Weibel, Phys. Rev. Lett. 2 83 (1959). 

[26] A. Stockem, M. E. Dieckmann, and R. Schlickeiser, 
Plasma Phys. Controll. Fusion 51 075014 (2009). 

[27] D. W. Forslund, and J. P. Freidberg, Phys. Rev. Lett. 27 
1189 (1971). 

[28] M. E. Dieckmann, H. Ahmed, G. Sarri, D. Doria, 1. 
Kourakis, L. Romagnani, M. Pohl, and M. Borghesi, 
Phys. Plasmas 20 042111 (2013). 

[29] J. F. Hansen, M. J. Edwards, D. H. Froula, A. D. Edens, 
G. Gregori, and T. Ditmire, Phys. Plasmas 13 112101 
(2006). 

[30] E. G. Berezhko, L. T. Ksenofontov, and H. J. Volk, As- 
tron. Astrophys. 412 Lll (2003). 

[31] A. R. Bell, Monthly Not. R. Astron. Soc. 353 550 (2004). 



